Esta Prueba de Evaluación Continuada cubre principalmente el módulo de generación de modelos no supervisados del programa de la asignatura.
En esta PEC trabajaremos la generación, interpretación y evaluación de un modelo de agregación k-means y otro DBSCAN. No perderemos de vista las fases de preparación de los datos, calidad del modelo y extracción inicial del conocimiento.
Material docente proporcionado por la UOC.
Módulo Métodos no supervisados del material didáctico.
Ejercicios teóricos
Todos los ejercicios deben ser presentados de forma razonada y clara, especificando todos y cada uno de los pasos que se hayan llevado a cabo para su resolución. No se aceptará ninguna respuesta que no esté claramente justificada.
Ejercicios prácticos
Para todas las PEC es necesario documentar en cada apartado del ejercicio práctico que se ha hecho y cómo se ha hecho.
El formato de entrega es: usernameestudiant-PECn.html (pdf o word) y rmd
Se debe entregar la PEC en el buzón de entregas del aula
A menudo es inevitable, al producir una obra multimedia, hacer uso de recursos creados por terceras personas. Es por lo tanto comprensible hacerlo en el marco de una práctica de los estudios de Informática, Multimedia y Telecomunicación de la UOC, siempre y cuando esto se documente claramente y no suponga plagio en la práctica.
Por lo tanto, al presentar una práctica que haga uso de recursos ajenos, se debe presentar junto con ella un documento en qué se detallen todos ellos, especificando el nombre de cada recurso, su autor, el lugar dónde se obtuvo y su estatus legal: si la obra está protegida por el copyright o se acoge a alguna otra licencia de uso (Creative Commons, licencia GNU, GPL …). El estudiante deberá asegurarse de que la licencia no impide específicamente su uso en el marco de la práctica. En caso de no encontrar la información correspondiente tendrá que asumir que la obra está protegida por copyright.
Deberéis, además, adjuntar los ficheros originales cuando las obras utilizadas sean digitales, y su código fuente si corresponde.
Los ejercicios se realizarán en base al juego de datos Hawks presente en el paquete R Stat2Data.
Los estudiantes y el profesorado del Cornell College en Mount Vernon, Iowa, recogieron datos durante muchos años en el mirador de halcones del lago MacBride, cerca de Iowa City, en el estado de Iowa. El conjunto de datos que analizamos aquí es un subconjunto del conjunto de datos original, utilizando sólo aquellas especies para las que había más de 10 observaciones. Los datos se recogieron en muestras aleatorias de tres especies diferentes de halcones: Colirrojo, Gavilán y Halcón de Cooper.
Hemos seleccionado este juego de datos por su parecido con el juego de datos penguins y por su potencial a la hora de aplicarle algoritmos de minería de datos no supervisados. Las variables numéricas en las que os basaréis son: Wing, Weight, Culmen, Hallux
library(reticulate)
use_condaenv("r-reticulate")Obtenemos el archivo .csv de datos de una fuente online alternativa, porque queremos trabajar con python y al pasar el datafame de R a un dataframe de Pandas en python la conversión de tipo de dato nos la hace mal.
import pandas as pd
raw_df = pd.read_csv('https://r-data.pmagunia.com/system/files/datasets/dataset-48225.csv')
raw_df.info()## <class 'pandas.core.frame.DataFrame'>
## RangeIndex: 908 entries, 0 to 907
## Data columns (total 19 columns):
## # Column Non-Null Count Dtype
## --- ------ -------------- -----
## 0 Month 908 non-null int64
## 1 Day 908 non-null int64
## 2 Year 908 non-null int64
## 3 CaptureTime 908 non-null object
## 4 ReleaseTime 907 non-null object
## 5 BandNumber 908 non-null object
## 6 Species 908 non-null object
## 7 Age 908 non-null object
## 8 Sex 332 non-null object
## 9 Wing 907 non-null float64
## 10 Weight 898 non-null float64
## 11 Culmen 901 non-null float64
## 12 Hallux 902 non-null float64
## 13 Tail 908 non-null int64
## 14 StandardTail 571 non-null float64
## 15 Tarsus 75 non-null float64
## 16 WingPitFat 77 non-null float64
## 17 KeelFat 567 non-null float64
## 18 Crop 565 non-null float64
## dtypes: float64(9), int64(4), object(6)
## memory usage: 134.9+ KBLas variables que vamos a tratar son las siguientes:
Wing: Longitud en mm de la pluma del ala primaria
Weight: Peso en gr
Culmen: Longitud en mm del pico superior
Hallux: Longitud en mm de la garra trasera
wanted_cols = ['Wing', 'Weight', 'Culmen', 'Hallux', 'Species']
data_cols = ['Wing', 'Weight', 'Culmen', 'Hallux']
df = raw_df[wanted_cols]
df.info()## <class 'pandas.core.frame.DataFrame'>
## RangeIndex: 908 entries, 0 to 907
## Data columns (total 5 columns):
## # Column Non-Null Count Dtype
## --- ------ -------------- -----
## 0 Wing 907 non-null float64
## 1 Weight 898 non-null float64
## 2 Culmen 901 non-null float64
## 3 Hallux 902 non-null float64
## 4 Species 908 non-null object
## dtypes: float64(4), object(1)
## memory usage: 35.6+ KBdf.dropna(inplace=True)## C:\Users\gevia\CONDA~1\envs\R-RETI~1\lib\site-packages\pandas\util\_decorators.py:311: SettingWithCopyWarning:
## A value is trying to be set on a copy of a slice from a DataFrame
##
## See the caveats in the documentation: https://pandas.pydata.org/pandas-docs/stable/user_guide/indexing.html#returning-a-view-versus-a-copy
## return func(*args, **kwargs)df.info()## <class 'pandas.core.frame.DataFrame'>
## Int64Index: 891 entries, 0 to 907
## Data columns (total 5 columns):
## # Column Non-Null Count Dtype
## --- ------ -------------- -----
## 0 Wing 891 non-null float64
## 1 Weight 891 non-null float64
## 2 Culmen 891 non-null float64
## 3 Hallux 891 non-null float64
## 4 Species 891 non-null object
## dtypes: float64(4), object(1)
## memory usage: 41.8+ KBdf.describe()## Wing Weight Culmen Hallux
## count 891.000000 891.000000 891.000000 891.000000
## mean 315.947475 771.615039 21.808923 26.413244
## std 95.316699 462.936901 7.293190 17.825541
## min 37.200000 56.000000 8.600000 9.500000
## 25% 202.000000 185.000000 12.800000 15.100000
## 50% 370.000000 970.000000 25.500000 29.400000
## 75% 390.000000 1120.000000 27.350000 31.400000
## max 480.000000 2030.000000 39.200000 341.400000Con los pasos anteriores nos hemos quedado sólo con las columnas de datos que nos interesan. La columna species nos puede servir de referencia en algún momento, pero no trabajaremos con ella porque vamos a tratar con algoritmos no supervisados, y por lo tanto no tenemos la etiqueta clasificatoria correcta. Simplemente trataremos el conjunto de datos y trataremos de agruparlos en distintos clusters. En este caso sabemos que el agrupamiento debería ser de 3, pero en un caso general no será así.
Además de quedarnos sólo con las columnas que nos interesan, también elminamos los registros con valores nulos, ya que son un número muy bajo y no merece la pena hacer otro tipo de tratamiento o imputación.
También podemos ver los valores descriptivos básicos. Para hacernos una mejor idea de la distribución de cada variable, realizamos un histograma.
import matplotlib.pyplot as plt
import seaborn as sns
sns.set_theme()
sns.set_palette('Set2')
fig=df.hist(bins=20)
plt.show()
Podemos observar que hay valores outliers en la variable Hallux. Eliminaremos esos registros de valores y nos quedaremos sólo con los que tienen un valor >50mm.
df.sort_values('Hallux').tail(10)## Wing Weight Culmen Hallux Species
## 360 478.0 1473.0 39.2 44.7 RT
## 265 387.0 1120.0 26.8 50.2 RT
## 740 260.0 565.0 19.7 54.5 CH
## 188 401.0 1405.0 29.1 82.8 RT
## 39 193.0 100.0 9.3 101.0 SS
## 582 170.0 110.0 11.4 121.0 SS
## 617 164.0 95.0 11.4 130.0 SS
## 102 198.0 188.0 12.2 143.0 SS
## 90 410.0 1210.0 27.5 308.0 RT
## 95 418.0 1180.0 30.1 341.4 RTdf = df[df['Hallux'] < 50]
df.describe()## Wing Weight Culmen Hallux
## count 882.000000 882.000000 882.000000 882.000000
## mean 316.222449 772.716553 21.830215 25.172676
## std 95.159444 462.097973 7.279885 8.222611
## min 37.200000 56.000000 8.600000 9.500000
## 25% 203.000000 185.000000 12.800000 15.025000
## 50% 370.000000 970.000000 25.500000 29.350000
## 75% 390.000000 1120.000000 27.300000 31.400000
## max 480.000000 2030.000000 39.200000 44.700000fig=df.hist(bins=20)
plt.show()
Ahora las distribuciones tienen más sentido y no se aprecian outliers graves.
Hacemos un gráfico de scatter con cada pareja de variables.
plt.figure()
fig = sns.pairplot(data=df, diag_kind='hist', hue='Species')
plt.show()
Vemos con los gráficos scatter de parejas de variables que en varios de los casos se distinguen tres clusters, si bien dos de ellos están bastante cerca. Además, podemos ver cómo todas las gráfica de pares de variables nos dan una información similar, por lo que todas ellas deben estar muy correladas. Lo comprobamos:
df.corr()## Wing Weight Culmen Hallux
## Wing 1.000000 0.934083 0.958302 0.934874
## Weight 0.934083 1.000000 0.952891 0.931660
## Culmen 0.958302 0.952891 1.000000 0.957087
## Hallux 0.934874 0.931660 0.957087 1.000000Efectivamente, la correlación es muy alta, especialmente entre las variables Wing, Weight y Cullmen.
Antes de realizar el clustering, vamos a comprobar si podemos reducir el número de variables del dataset para que la representación visual de los datos sea mucho más sencilla. Para ello, vamos a hacer uso de la descomposición PCA (principal components analysis) y comprobaremos sin con sólo dos de las componentes resultantes perdemos o no mucha de la información original.
from sklearn.decomposition import PCA
from sklearn.preprocessing import StandardScaler
X = StandardScaler().fit_transform(df[data_cols])
n_pca = 4
pca = PCA(n_components=n_pca)
pca_val = pca.fit_transform(X)
pca_df = pd.DataFrame(pca_val, columns=['pca1', 'pca2', 'pca3', 'pca4'])
# Unimos al DataFrame con las componentes PCA la etiqueta de especie como referencia
final_df = pca_df.assign(Species=df['Species'].values)
print(final_df.head())## pca1 pca2 pca3 pca4 Species
## 0 1.086921 0.224645 -0.182713 -0.042647 RT
## 1 1.283533 0.201659 -0.030918 0.009588 RT
## 2 -0.913645 0.305847 0.132601 0.018886 CH
## 3 -2.539926 0.005589 -0.119992 -0.006538 SS
## 4 1.733411 0.147379 -0.174137 0.031517 RTpca_var = pca.explained_variance_ratio_
plt.figure()
pfig = plt.bar(x=range(n_pca), height=pca_var)
plt.xlabel('PCA component')
plt.ylabel('Explained variance ratio')
plt.xticks(list(range(n_pca)))## ([<matplotlib.axis.XTick object at 0x000000004F02A760>, <matplotlib.axis.XTick object at 0x000000004F02A5E0>, <matplotlib.axis.XTick object at 0x000000004F04D910>, <matplotlib.axis.XTick object at 0x00000000500B09D0>], [Text(0, 0, ''), Text(0, 0, ''), Text(0, 0, ''), Text(0, 0, '')])for i, v in enumerate(pca_var):
plt.gca().text(i-0.09, v+0.01, f'{v:.2f}', fontweight='bold')
plt.show()
Vemos en la gráfica cómo sólo la primera de las componentes PCA contiene el 96% de la información original. Es un porcentaje tan alto que nos podríamos quedar sólo con esa dimensión, pero tomaremos las dos primeras al tener espacio para representar 2 en 2D.
Dibujamos en el mapa 2D las muestras que tenemos según sus dos primeras componentes principales. Añadimos en color la información de la clase a la que sabemos que pertenece.
plt.figure()
sns.scatterplot(x='pca1', y='pca2', hue='Species', data=final_df)
plt.title('Hawks samples by species')
plt.show()
De esta respresentación se pueden extraer varias conclusiones:
Efectivamente, con una sóla componente PCA seríamos capaces de hacer una separación muy buena entre clases.
Visualmente se puede ver cómo las muestras no están perfectamente separadas.
Eliminando la información de la clase (color), se puede ver cómo hay un claro cluster en la derecha y a la izquierda tenemos 2 o 4, según seamos de estrictos con su definición.
Hay tres clases, por lo que cabría esperar una agrupación en tres clusters, siempre que dispongamos de esta información de antemano, que en muchos casos no es así.
No todos los clusters que se aprecian visualmente tienen el mismo tamaño, el mismo número de muestras, ni la misma densidad.
Vamos a aplicar ahora el algoritmo KMeans sobre esta versión reducida a dos componentes PCA, pues sabemos que los resultados serán muy similares a si utilizásemos las cuatro otriginales. Hacemos un barrido con distinto número de clusters: entre 2 y 8.
from sklearn.cluster import KMeans
distortions = []
for n_clusters in range(2,9):
k_means = KMeans(n_clusters=n_clusters, random_state=41)
k_means.fit(final_df[['pca1', 'pca2']])
distortions.append(k_means.inertia_)## KMeans(n_clusters=2, random_state=41)
## KMeans(n_clusters=3, random_state=41)
## KMeans(n_clusters=4, random_state=41)
## KMeans(n_clusters=5, random_state=41)
## KMeans(n_clusters=6, random_state=41)
## KMeans(n_clusters=7, random_state=41)
## KMeans(random_state=41)plt.figure()
fig=plt.plot(range(2,9), distortions)
plt.xlabel('Number of clusters')
plt.ylabel('SSE')
plt.show()
En la gráfica se muestra el error cuadrático medio (SSE) de las muestras dentro de su cluster. A menor valor, mejor es la agrupación en clusters. Como se puede observar, la curva es descendiente, ya que el error siempre va a disminuir al aumentar el número de clusters. Aplicando el método elbow para determinar el número adecuado de clusters, vemos que el valor óptimo de clusters estaría en 4, ya que de 3 a 4 el SSE disminuye mucho, pero a partir de ese punto el SSE se estabiliza bastante y disminuye ya de forma mucho más lenta. Vemos, sin embargo, que este número no coincide con el número de clases de halcones que sabemos que existen en el dataset de antemano, por lo que con este algorito y este método de evaluación la elección del número clusters sin información previa no sería correcta.
Vamos a representar gráficamente la asignación de las muestras en cuatro clusters, el número que hemos obtenido como óptimo con el método anterior. Dibujamos también los centroides de cada uno de los clusters.
n_clusters = 4
k_means = KMeans(n_clusters=n_clusters, random_state=41)
classes = k_means.fit_predict(final_df[['pca1', 'pca2']])
plt.figure()
sns.scatterplot(x='pca1', y='pca2', data=final_df, hue=classes, palette='Set2')
for i in range(n_clusters):
plt.plot(
k_means.cluster_centers_[i, 0],
k_means.cluster_centers_[i, 1],
marker='o',
markerfacecolor='r',
markeredgecolor='k',
markersize=8)
plt.show()
Vemos cómo el gran cluster de la derecha se ha dividido en dos, cuando la lógica nos llevaría a pensar que debería ser uno sólo. Esto implica que el método utilizado no nos ha dado un resultado satisfactorio.
Vamos a probar con 3 clusters, el número de clases que sabemos que existen de antemano.
n_clusters = 3
k_means = KMeans(n_clusters=n_clusters, random_state=41)
classes = k_means.fit_predict(final_df[['pca1', 'pca2']])
plt.figure()
sns.scatterplot(x='pca1', y='pca2', data=final_df, hue=classes, palette='Set2')
for i in range(n_clusters):
plt.plot(
k_means.cluster_centers_[i, 0],
k_means.cluster_centers_[i, 1],
marker='o',
markerfacecolor='r',
markeredgecolor='k',
markersize=8)
plt.show()
Se puede ver ahora cómo la agrupación en los 3 clusters que veíamos visualmente es muy buena, por lo que este valor hubiese sido una elección mucho más adecuada. Hay que tener en cuenta que el agrupamiento realizado no nos da como resultado la clase a la que pertenece la muestra (en este caso, especie de halcón), sino sólamente diferencia entre tres clases distintas.
Comparamos ahora esta última imagen con las clases reales que teneníamos y comprobamos, que, efectivamente, esta agrupación habría sido mucho mejor y habría coincidido mucho más con la agrupacion que hemos realizado con tres clusters. Evidentemente, las muestras que estaban muy lejos de su cluster, serían mal agrupadas.
plt.figure()
fig = sns.scatterplot(x='pca1', y='pca2', data=final_df, hue= 'Species')
for i in range(n_clusters):
plt.plot(
k_means.cluster_centers_[i, 0],
k_means.cluster_centers_[i, 1],
marker='o',
markerfacecolor='r',
markeredgecolor='k',
markersize=8)
plt.show()
Realizamos el análists utilizando una métrica diferente, con el coeficiente silhoutte. Esta métrica tiene en cuenta no sólamente la cercanía a su centroide, sino también la lejanía al resto de centroides de otros clusters, de forma que la primera pueda ser minimizada y la segunda maximizada. Su valor estará siempre entre -1 y 1, donde 1 será una clasificación óptima en la muestra está justo en el centroide de su cluster y muy separada del resto de centroides. El coeficiente silhouette se calcula para cada una de las muestras, a partir del cual se puede obtener la media para todas las muestras del dataset.
import numpy as np
from sklearn.metrics import silhouette_samples, silhouette_score
range_n_clusters = [2, 3, 4, 5]
for n_clusters in range_n_clusters:
fig, (ax1, ax2) = plt.subplots(1, 2)
fig.set_size_inches(22, 10)
fig.suptitle(f'n_clusters = {n_clusters}', fontsize=14)
ax1.set_xlim([-0.1, 1])
ax1.set_ylim([0, final_df.shape[0] + (n_clusters + 1) * 10])
clusterer = KMeans(n_clusters=n_clusters, random_state=10)
cluster_labels = clusterer.fit_predict(final_df[['pca1', 'pca2']])
silhouette_avg = silhouette_score(final_df[['pca1', 'pca2']], cluster_labels)
sample_silhouette_values = silhouette_samples(final_df[['pca1', 'pca2']], cluster_labels)
y_lower = 10
for i in range(n_clusters):
ith_cluster_silhouette_values = sample_silhouette_values[cluster_labels == i]
ith_cluster_silhouette_values.sort()
size_cluster_i = ith_cluster_silhouette_values.shape[0]
y_upper = y_lower + size_cluster_i
ax1.fill_betweenx(np.arange(y_lower, y_upper), 0, ith_cluster_silhouette_values, alpha=0.7)
ax1.text(-0.05, y_lower + 0.5 * size_cluster_i, str(i))
y_lower = y_upper + 10 # 10 for the 0 samples
ax1.set_title('Silhouette plot for the clusters')
ax1.set_xlabel('The silhouette coefficient')
ax1.set_ylabel('Cluster label')
ax1.axvline(x=silhouette_avg, color='red', linestyle='--')
ax1.set_xticks([0, 0.2, 0.4, 0.6, 0.8, 1])
sns.scatterplot(x='pca1', y='pca2', data=final_df, ax=ax2, hue=cluster_labels, palette='Set2')
centers = clusterer.cluster_centers_
# Draw white circles at cluster centers
ax2.scatter(
centers[:, 0],
centers[:, 1],
marker="o",
c="white",
alpha=1,
s=200,
edgecolor="k",
)
for i, c in enumerate(centers):
ax2.scatter(c[0], c[1], marker='$%d$' % i, alpha=1, s=50, edgecolor="k")
ax2.set_title("Clustered data")
ax2.set_xlabel("PCA1")
ax2.set_ylabel("PCA2")
plt.show()
(para esta última parte se ha utilizado como referencia el siguiente análisis : https://scikit-learn.org/stable/auto_examples/cluster/plot_kmeans_silhouette_analysis.html)
En este caso lo que hecemos es dibujar, para cada valor de n_clusters, el coeficiente de silhoutte de todas las muestras por cluster junto con su asignación por muestras en el espacio 2D en base a sus dos componentes principales PCA. Por otro lado, en las gráficas de las izquierda se dibuja también una gráfica vertical roja que corresponde a la media del coeficiente silhoutte de todas las muestras. Este número nos da bastante información y puede tomarse como métrica como hemos hecho anteriormente con el error cuadrático medio, pero pormenorizar muestra a muestra y dividir por clusters puede darnos mucha más información. Analizando los gráficos anteriores podemos extraer lo siguiente:
Con dos clusters, la separación es la que podríamos esperar a partir del scatter plot y la media del coeficiente silhouette es bastante alta, favorecida por el hecho de que casi todas las muestras de ambos clusters se encuentran muy separadas del centroide del cluster que no es el suyo. Podemos observar en el perfil del coeficiente silhouette de la clase 1 (naranja) como hay un buen número que bajan mucho respecto a la media, que correponden con las muestras de la derecha que se encuentran mucho más cerca del otro cluster (verde). La mayoría de muestras de esta clase 1 tiene un coeficiente muy alto porque en realidad la zona cercana al centroide de dicho cluster es mucho más densa.
Con tres clusters, tenemos la agrupación que nos esperaríamos, pero vemos cómo uno de los clusters, el 2 (azul), además de tener muy pocas muestras en comparación con los otros dos, tiene todas sus muestras con un coeficiente silhouette muy por debajo de la media. Esto suele ser un signo de que dicho cluster debería sobrar porque está o bien muy disperso o bien muy cerca a otro cluster.
Con cuatro clusters, el perfil de los coeficientes de las tres clases tiene mejor aspecto porque todos superan la media, pero vemos que ésta ha bajado mucho (<0.6) respecto al caso de dos clusters (>0.8), por lo que este último sería un valor más adecuado para escoger. Además, podemos observar en el perfil de las clases 1 y 3 que tienen una pendiente mucho menos pronunciada, lo que quiere decir que estas clases están mucho más dispersas. Si vemos las distribución en scatter, vemos cómo estos dos clusters se corresponden con el gran cluster de la derecha, que se ha dividido en dos.
Con cinco clusters, tenemos el mismo problema que con cuatro, y un coeficiente de silhoutte medio muy similar. En este caso, el cluster de más a la izquierda se ha dividido en dos, si bien esta división sí que parece bastante acertada a juzgar por los perfiles del coeficiente y la representación en scatter de los clusters.
En base a todo lo anterior, con esta métrica elegiríamos un número n_clusters=2, a diferencia de lo que haríamos con la evaluación de SSE.
Vamos a continuar el ejercicio 2 a partir de la transformación PCA que hemos realizado en el ejercicio anterior teniendo en cuenta que con las dos primeras componentes tenemos un altísimo 98% de la información y de esta manera podemos representar de mejor forma los puntos en el espacio 2D.
Aplicamos en primer lugar el algoritmo DBSCAN. Para ello, debemos especificar los parámetros eps y min_samples.
from sklearn.cluster import DBSCAN
eps = 0.3
min_samples = 10
db = DBSCAN(eps=eps, min_samples=min_samples).fit(final_df[['pca1', 'pca2']])
labels = db.labels_
if -1 in labels:
n_clusters = len(set(labels)) - 1
else:
n_clusters = len(set(labels))
print(f'Number of clusters: {n_clusters}')## Number of clusters: 3plt.figure()
palette ={-1: "k", 0:'b', 1:'r', 2:'g', 3:'gold', 4:'violet', 5:'teal'}
fig = sns.scatterplot(x='pca1', y='pca2', data=final_df, hue=labels, palette=palette)
plt.show()
Vemos que el algoritmo en este caso nos devuelve 3 clusters, y este número no ha sido especificado de antemano como en el caso de kmeans, sino que es obtenido a partir de los parámetros eps y min_samples especificados. Además, nos devuelve una cuarta clase (etiquetada como -1, en color negro) que contiene los elementos que ha identificado como outliers. Vemos también que la agregación en estos tres clusters es bastante buena, ya que distribuye los clusters como lo podríamos hacer de forma visual.
A continuación, probamos otros valores de eps.
eps = 0.5
min_samples = 10
db = DBSCAN(eps=eps, min_samples=min_samples).fit(final_df[['pca1', 'pca2']])
labels = db.labels_
if -1 in labels:
n_clusters = len(set(labels)) - 1
else:
n_clusters = len(set(labels))
print(f'Number of clusters: {n_clusters}')## Number of clusters: 2plt.figure()
fig = sns.scatterplot(x='pca1', y='pca2', data=final_df, hue=labels, palette=palette)
plt.show()
eps = 0.2
min_samples = 10
db = DBSCAN(eps=eps, min_samples=min_samples).fit(final_df[['pca1', 'pca2']])
labels = db.labels_
if -1 in labels:
n_clusters = len(set(labels)) - 1
else:
n_clusters = len(set(labels))
print(f'Number of clusters: {n_clusters}')## Number of clusters: 4plt.figure()
fig = sns.scatterplot(x='pca1', y='pca2', data=final_df, hue=labels, palette=palette)
plt.show()
De las gráficas anteriores vemos cómo si aumentamos el valor de eps obtenemos un número menor de clusters, uniendo los dos que quedaban a la izquierda. Por otro lado, disminuyendo su valor aumentamos el número de clusters a cuatro, separando en dos clusters uno de los iniciales. El cambio en este valor también afecta a los outliers, pues un eps menor detectará más elementos como outliers al tener los clusters de menor tamaño, como se puede comprobar en la última gráfica.
import numpy as np
from sklearn.metrics import calinski_harabasz_score
silhouette = []
ch_score = []
for eps_val in np.linspace(0.1,0.8,8):
db = DBSCAN(eps=eps_val, min_samples=min_samples).fit(final_df[['pca1', 'pca2']])
ch_score.append(calinski_harabasz_score(final_df[['pca1', 'pca2']], db.labels_))
silhouette.append(silhouette_score(final_df[['pca1', 'pca2']], db.labels_))
plt.figure()
fig=plt.plot(np.linspace(0.1,0.8,8), silhouette)
plt.xlabel('eps')
plt.ylabel('Silhouette score')
plt.show()
Se puede ver cómo el coeficiente de silhouette promedio mejora con el incremento de la variable eps. Si aplicamos el método elbow sobre esta gráfica, elegiríamos un valor eps de 0.3, ya que a partir de ese punto la mejora no es muy notable. Este valor de 0.3 coincide justo con el valor que hemos elegido en el primer caso, y que como hemos visto acaba creando tres clusters que agregan las muestras bastante bien.
Hacemos el mismo estudio más detallado que hemos hecho con el coeficiente silhouette en el caso de kmeans, pero con DBSCAN, haciendo esta vez barrido en eps en lugar de n_clusters.
import numpy as np
from sklearn.metrics import silhouette_samples, silhouette_score
eps_list = list(np.linspace(0.1,0.5,5))
for eps in eps_list:
fig, (ax1, ax2) = plt.subplots(1, 2)
fig.set_size_inches(18, 10)
fig.suptitle(f'eps = {eps:.1f}', fontsize=14)
ax1.set_xlim([-0.1, 1])
ax1.set_ylim([0, final_df.shape[0] + (n_clusters + 1) * 10])
clusterer = DBSCAN(eps=eps, min_samples=min_samples).fit(final_df[['pca1', 'pca2']])
cluster_labels = clusterer.fit_predict(final_df[['pca1', 'pca2']])
silhouette_avg = silhouette_score(final_df[['pca1', 'pca2']], cluster_labels)
sample_silhouette_values = silhouette_samples(final_df[['pca1', 'pca2']], cluster_labels)
y_lower = 10
if -1 in labels:
n_clusters = len(set(labels)) - 1
else:
n_clusters = len(set(labels))
for i in range(n_clusters):
ith_cluster_silhouette_values = sample_silhouette_values[cluster_labels == i]
ith_cluster_silhouette_values.sort()
size_cluster_i = ith_cluster_silhouette_values.shape[0]
y_upper = y_lower + size_cluster_i
ax1.fill_betweenx(np.arange(y_lower, y_upper), 0, ith_cluster_silhouette_values, alpha=0.7, color=palette[i])
ax1.text(-0.05, y_lower + 0.5 * size_cluster_i, str(i))
y_lower = y_upper + 10 # 10 for the 0 samples
ax1.set_title('Silhouette plot for the clusters')
ax1.set_xlabel('The silhouette coefficient')
ax1.set_ylabel('Cluster label')
ax1.axvline(x=silhouette_avg, color='red', linestyle='--')
ax1.set_xticks([0, 0.2, 0.4, 0.6, 0.8, 1])
sns.scatterplot(x='pca1', y='pca2', data=final_df, ax=ax2, hue=cluster_labels, palette=palette)
ax2.set_title("Clustered data")
ax2.set_xlabel("PCA1")
ax2.set_ylabel("PCA2")
plt.show()
Podemos sacar las siguientes conclusiones:
Con eps=0.1 tenemos demasiadas muestras cataligadas como outliers, además de un cluster demasiado pequeño.
Con eps=0.2 mejoramos bastante la media del coeficiente, pero tenemos dos clusters demasiado pequeños.
Con eps=0.3 tenemos un buen nivel de coeficiente silhoutte medio, además de unos buenos perfiles de clusters, todos ellos superando ese nivel medio y con una cantidad de muestras adecuada.
Con eps=0.4 y más, obtenemos sólo dos agrupaciones y el coeficiente medio es también bueno, pero se puede apreciar cómo en la clase 1 (roja), hay un buen número de muestras con un coeficiente muy bajo, por lo que la elección de eps=0.3 parece más acertada.
En base a lo anterior, y en consonancia con lo visto en la gráfica de coeficiente silhouette medio vs eps, elegiríamos un valor de eps=0.3.
Utilizamos a continuación un índice distinto para la evaluación del modelo, en este caso el Calinski and Harabasz.
plt.figure()
fig=plt.plot(np.linspace(0.1,0.8,8), ch_score)
plt.xlabel('eps')
plt.ylabel('Calisnki-Harabasz score')
plt.show()
Aun tratándose de un índice distinto, vemos que el resultado que obtenemos es el mismo e incluso de forma más clara que con el índice silhouette: el valor de eps óptimo siguiendo el método elbow es de 0.3.
Ahora vamos a aplicar el algoritmo OPTICS con los mismos datos. El algoritmo OPTICS funciona de forma muy similar al DBSCAN, pero en este caso no tenemos que fijar el valor de eps de antemano.
from sklearn.cluster import OPTICS
min_samples = 60
clust = OPTICS(min_samples=min_samples).fit(final_df[['pca1', 'pca2']])
labels = clust.labels_
if -1 in labels:
n_clusters = len(set(labels)) - 1
else:
n_clusters = len(set(labels))
print(f'Number of clusters: {n_clusters}')## Number of clusters: 3plt.figure()
fig = sns.scatterplot(x='pca1', y='pca2', data=final_df, hue=labels, palette=palette)
plt.show()
Vemos como con el mismo set de datos este algoritmo nos revuelve un agrupamiento muy distinto al de DBSCAN, aun resultando en el mismo número de clusters. Esto se debe a que el valor de min_samples especificado es demasiado alto. Vamos a ver lo que ocurre cuando lo disminuimos:
min_samples = 30
clust = OPTICS(min_samples=min_samples).fit(final_df[['pca1', 'pca2']])
labels = clust.labels_
if -1 in labels:
n_clusters = len(set(labels)) - 1
else:
n_clusters = len(set(labels))
print(f'Number of clusters: {n_clusters}')## Number of clusters: 4plt.figure()
fig = sns.scatterplot(x='pca1', y='pca2', data=final_df, hue=labels, palette=palette)
plt.show()
Vemos ahora cómo ahora ha aparecido un cluster más, pero la agrupación sigue sin ser la que esperaríamos, y hay un gran número de muestras catalogadas como outliers.
EL algoritmo OPTICS en realidad asigna a cada muestra una reachability distance, que será la mayor valor entre la core distance y la distancia euclídea entre la muestra y su core. Dibujamos esta reachability distance ordenando los puntos y asignando un color a cada clase, donde los puntos negros son outliers.
reachability = clust.reachability_[clust.ordering_]
labels = clust.labels_[clust.ordering_]
space = np.arange(len(labels))
plt.figure()
for cl in range(0, 5):
X = space[labels == cl]
Y = reachability[labels == cl]
plt.plot(X, Y, color=palette[cl], alpha=0.3, linestyle='', marker='o', markersize=3)
plt.plot(space[labels == -1], reachability[labels == -1], "k.", alpha=0.3)
#plt.plot(space, np.full_like(space, 0.5, dtype=float), "k-.", alpha=0.5)
plt.ylabel("Reachability")
plt.show()
Vemos como en el reachability plot existen zonas valle que definen cada uno de los clusters con su color correspondiente. En negro se mantienen las muestras clasificadas como outliers. Además, cuanto más profundo el valle querrá decir que mayor es la densidad de puntos de ese cluster. Esto lo podemos comprobar porque los clusters verde y amarillo son los más densos y son en los que el valle es más profundo; el caso opuesto lo tenemos con el cluster rojo, que tiene el valle más alto de todos, y se puede comprobar claramente en el scatter plot de muestras que es el que tiene las muestras más dispersas.
Hemos realizado un análisis y comparativa de los algoritmos kmeans y DBSCAN para el clustering de datos. Con el algoritmo kmeans, si partimos de la información previa de que existen tres clases, y por lo tanto fijamos n_clusters=3, podemos obtener un agrupamiento muy bueno en tres clusters de acuerdo con lo que visualmente podemos apreciar y con las etiquetas con las que partimos. Sin embargo, si no partimos de esta información, la elección sería la de cuatro clusters si utilizamos la métrica de error cuadrático medio o de dos si utilizamos el coeficiente silhouette, en ningún caso el número de clases que sabemos que existen.
Con el algoritmo DBSCAN, sin embargo, obtenemos un buen agrupamiento con en tres clusters sin disponer de esta información de partida, y obtenemos el mismo resultado con las dos métricas de evaluación que hemos probado (coeficiente sihouette y Calinski-Harabasz) . Con este algoritmo, eso sí, hay que tener especial cuidado con la elección de los parámetos que lo definen porque de ello depende que tengamos un buen clustering o no. Además, tenemos la ventaja adicional de poder detectar datos outlier.
Por lo tanto, con este conjunto de datos nos quedaríamos con el algoritmo DBSCAN como el mejor para realizar un clustering adecuado.
Las ventajas de kmeans son las siguientes:
Es mucho más rápido, lo que permite utilizarlo cuando el número de puntos es grande.
Sólo requiere definir un parámetro.
No se ve afectado por variaciones en la densidad de los clusters.
Los inconvenientes de kmeans son los siguientes:
Los clusters tienen una forma aproximadamente esférica, por lo que no es válido si los clusters no tienen dicha forma.
Tienes que definir de antemano el número de clusters a utilizar.
Es muy sensible a la presencia de outliers en los datos.
Las ventajas de DBSCAN son las siguientes:
Detecta clusters de forma arbitraria, a diferencia de k-means.
No necesita que se le especifique el número de clusters a priori, algo que k-means sí necesita.
Es menos sensible a ruido y robusto frente a la presencia de valores outliers.
Permite la detección directa de valores outliers dentro de los datos.
Los inconvenientes de DBSCAN son:
La elección de los parámetros eps y min_samples puede ser compleja, como hemos comprobado en el ejercicio.
No consigue realizar un buen agrupamiento cuando hay clusters con densidades de puntos muy distintas.
En el trabajo se propone la agregación en clusters de muestras de medidas de cuatro parámetros de varias especies de halcones. En primer lugar, se han eliminado muestras con valores nulos y eliminado algunos outliers claros. A continuación, se ha aplicado la técnica PCA de reducción de dimensionalidad para poder representar las muestras en un espacio 2D y comprobado que con ello no hemos perdido casi información.
A continuación se ha podido comprobar el funcionamiento del algoritmo k-means con distinto número de clusters y evaluado con distintas métricas. Hemos realizado también una comparativa con las etiquetas que teníamos de las especies previamente, y hemos comprobado como escogiendo el número de clusters igual al número de clases, podemos hacer un agrupamiento que coincide muy bien con el de la clasificación por especies.
Posteriormente, hemos realizado un estudio similar aplicando el algoritmo DBSCAN y esta vez viendo como cambia su funcionamiento variando el parámetro eps y evaluando el modelo con distintas métricas. Se interpreta los clusters obtenidos en los distintos casos, incluyendo los outliers